*! version 1.0  27feb2020

version 10.1

program define warofatt, eclass

  syntax varlist [if] [in], [TSTarvalues(numlist) TSTARMIN(numlist max=1) TSTARMAX(numlist max=1) PRobss(varlist) ///
	iterate(numlist int max=1) SAVEloglik(string) PLOTLoglikelihood GRSAVELoglikelihood(string) MLmaxoptions(string) ///
	PLOTHazard GRSAVEHazard(string) STartvalues ROBust CLuster(varname) Level(int $S_level) *  ]

		/* This is probably not the best way to pass these, but I can't find an equivalent 	*/
		/* to -mlopts- for noninteractive mode options like init() and search(off). This works, */
		/* but will allow unacceptable options to be passed to the optional -twoway- commands.	*/
	
	mlopts mlopts rest, `options'

	if "`iterate'" != "" {
	  local maxiters iterate(`iterate')
	  }

	else if "`iterate'" == "" {
	  local iterate = c(maxiter)
	  }
	  
	_get_gropts, graphopts(`rest') gettwoway
        local gropts `"`s(twowayopts)'"'
	
	marksample touse
	
	if "`s(graphopts)'" != "" {
	  di in red _n "warofatt for Stata does not yet support the `s(graphopts)' option"
	  exit 198
	  }
	
	       /* Estimate robust SEs. */
	
	
	if "`robust'" != "" { 
		local robust robust 
		}
		
	       /* Estimate clustered SEs. */
	
	if "`cluster'" != "" { 
	  local clusopt "cl(`cluster')"
	  if "`robust'" != "" {
		di in red "Error: the robust and cluster options are incompatible"
		exit 198
		}
	  }
	  
	if "`iterate'" != "" {
	  local maxiters iterate(`iterate')
	  }

	else if "`iterate'" == "" {
	  local iterate = c(maxiter)
	  }
	  
	gettoken depvar vars_W : varlist
	  
	  
		/****************************************************/
		/* Get values of t* to search over if not provided. */
		/****************************************************/


	if "`tstarvalues'" == "" {
	
	  quietly levelsof `depvar', local(tstars)
	  
	  quietly summarize `depvar'
	  
	  if "`tstarmin'" == "" {  
	  
		local tstarmin = r(min)
	  
		}
	  
	  if "`tstarmax'" == "" {  
	  
		local tstarmax = r(max)
		
		}
	  
	  }
	  
	else {

	  if "`tstarmin'" != "" {  
	  
		display in red "Error: can not specify both tstarvalues and tstarmin options."
		exit 198
	  
		}
	  
	  else if "`tstarmax'" != "" {  
	  
		display in red "Error: can not specify both tstarvalues and tstarmax options."
		exit 198
		
		}
	  
	  else {
	  
		local tstars `tstarvalues'
		local tstarmin = 0
		local tstarmax = .
		
		}
	  
	  }
	  
	local tstars_num = wordcount("`tstars'")

		/* Create variables to store likelihoods. */
	  
	tempvar tstar_var ll_tstar ll_value ll_convg ll_iters
	  
	quietly generat `ll_tstar' = .
	quietly generat `ll_value' = .
	quietly generat `ll_convg' = .
	quietly generat `ll_iters' = .

	quietly generat `tstar_var' = .


		/*********************************************/
		/* Run the model for specified values of t*. */
		/*********************************************/
	
	
	display ""
	noisily _dots 0, title("Performing MLE for requested values of t*")
	
	local i = 1
  
	foreach tstar of local tstars {
	
	  if `tstar' >= `tstarmin' & `tstar' <= `tstarmax' {

		  quietly replace `tstar_var' = `tstar'

		  if "`startvalues'"!="" { 
		  
			if regexm("`mlmaxoptions'","init") {
			  di in red "Error: can not specify both startvalues and ml maximize's init() option."
			  exit 198 
			  }

			warofatt_sv, tstar(`tstar') depvar(`depvar') indvars(`vars_W') probss(`probss') 
			  matrix _b0=r(b0)
			  local stvals = "init(_b0, copy) search(off)"
			  
			}
			  
		  capture ml model lf warofatt_llf (Weak: `depvar' `tstar_var' = `vars_W') (StrongStrong: ) ///
			(logit_pi: `probss') /ln_p_W  /ln_p_SS if `touse', ///
			maximize `robust' `clusopt' nooutput `maxiters' `stvals' `mlopts' `mlmaxoptions'

			if _rc == 198 {
			
			   di in red "Error: one of the following options not allowed by maximize:"
			   di in red "     `mlmaxoptions'"
			   exit 198
			
			  }
			  
			else if _rc == 0 {
			
			  quietly replace `ll_tstar' = `tstar' if _n == `i'
			  quietly replace `ll_value' = e(ll) if _n == `i'
			  quietly replace `ll_convg' = e(converged) if _n == `i'
			  quietly replace `ll_iters' = e(ic) if _n == `i'

			  noisily _dots `i' 0
			
			  }
			  
			else {
			
			  quietly replace `ll_tstar' = `tstar' if _n == `i'
			  quietly replace `ll_value' = . if _n == `i'
			  quietly replace `ll_convg' = . if _n == `i'
			  quietly replace `ll_iters' = . if _n == `i'

			  noisily _dots `i' 1

			  }
			  
			local i = `i'+1
			
		  }
		  
		}


		/***************************************************/
		/* Output the log-likelihood results if requested. */
		/***************************************************/

		
	if "`saveloglik'" != "" {
	
	  preserve
	  
	  quietly keep `ll_tstar' `ll_value' `ll_convg' `ll_iters'
	  quietly drop if missing(`ll_tstar')
	  
	    rename `ll_tstar' ll_tstar
	    rename `ll_value' ll_value
	    rename `ll_convg' ll_convg
	    rename `ll_iters' ll_iters
	    generat ll_maxiters = `iterate'
	  
	  quietly save `saveloglik'
	  
	  restore
	  
	  }

	  
		/**********************************************************/
		/* Rerun for the value that maximizes the log-likelihood. */
		/**********************************************************/

		
	display ""
	
	quietly sum `ll_value' if `ll_convg' == 1 & `ll_iters' < `iterate'

	  quietly sum `ll_tstar' if `ll_value' == r(max)
	
		local tstar_max = r(min)
		
		quietly replace `tstar_var' = `tstar_max'

		  if "`startvalues'"!="" {
		  
			if regexm("`mlmaxopts'","init") {
			  di in red "Error: can not specify both startvalues and ml maximize's init() option."
			  exit 198 
			  }
			  
			warofatt_sv, tstar(`tstar_max') depvar(`depvar') indvars(`vars_W') probss(`probss') 
			  matrix _b0=r(b0)
			  local stvals = "init(_b0, copy) search(off)"
			  
			}

		  ereturn clear
		  
		  ml model lf warofatt_llf (Weak: `depvar' `tstar_var' = `vars_W') (StrongStrong: ) (logit_pi: `probss') /ln_p_W /ln_p_SS if `touse', ///
			maximize `robust' `clusopt' nooutput iterate(`iterate') `stvals' `mlopts' `mlmaxoptions' ///
			title( "War of attrition duration model, t*=`tstar_max'")
			
		local F_W_tstar = 1 - exp(-(exp([Weak]_b[_cons])*`tstar_max')^exp([ln_p_W]_b[_cons]))

		di in r "" 
		di in r "" 
		
		ml display, level(`level') neq(3) plus
		_diparm ln_p_W, level(`level') prob label("ln_p_W")
		_diparm ln_p_W, level(`level') prob label("p_W") exp
		_diparm __bot__
		_diparm ln_p_SS, level(`level') prob label("ln_p_SS")
		_diparm ln_p_SS, level(`level') prob label("p_SS") exp
		_diparm __bot__
		di in smcl in gr "(Nontruncated) Weibull density evaluated at t*: " `F_W_tstar'
		di in smcl in gr "{hline 78}"

				
		/*********************************************/
		/* Graph the likelihood maxima if requested. */
		/*********************************************/

		
	if "`plothazard'" != "" & "`plotloglikelihood'" != "" {
	
	  di in red "Error: may not specify both the plotloglikelihood and plothazard options."
	  exit 198 
	
	  }

	
	else if "`plotloglikelihood'" != "" {

	  twoway connected `ll_value' `ll_tstar' if `ll_convg'==1 & `ll_iters' < `iterate', ///
		xtitle(Value of t*) ///
		ytitle(Final Log-likelihood) ///
		`gropts' ///
		saving(`grsaveloglikelihood')
	  
	  }


		/*********************************************/
		/* Graph the estimated hazard if requested. */
		/*********************************************/

	
	else if "`plothazard'" != "" {
		
	  quietly summarize `depvar' if e(sample)
	  
		local depvar_max=r(max)
		
	  local prcure : display %3.2f (1+exp(-[logit_pi]_b[_cons]))^(-1)
	  local ll : display %3.2f e(ll)
	  local tstar_max : display %3.2f `tstar_max'
  
	  local lambda_W = exp([Weak]_b[_cons])
	  local p_W = exp([ln_p_W]_b[_cons])
		
	  local lambda_SS = exp([StrongStrong]_b[_cons])
	  local p_SS = exp([ln_p_SS]_b[_cons])

		
	  twoway function `lambda_SS'*`p_SS'*((`lambda_SS'*x)^(`p_SS'-1)), range(0 `depvar_max') scheme(s1mono) ///
		lcolor(gs0) lwidth(thick) lpattern(solid) ///
	    || 	function y=(`lambda_W'*`p_W'*(`lambda_W'*x)^(`p_W'-1))*((exp(-(`lambda_W'*x)^`p_W'))/(exp(-(`lambda_W'*x)^`p_W') - exp(-(`lambda_W'*`tstar_max')^`p_W'))), ///
		range(0 `=`tstar_max' - 1') ///
		lcolor(gs0) lwidth(thick) lpattern(solid) ///
		title( "t*=`tstar_max' Pr(Both Players Strong)= `prcure' ll=`ll'") ///
		ylabel(#5, grid) ytitle("") ///
		xtitle("") ///
		xline(`tstar_max', lcolor(gs10) lwidth(medthin)) ///
		legend(off) ///
		`gropts' ///
		saving(`grsavehazard')
		
	  }
		
	ereturn scalar tstarmax = `tstar_max'
	ereturn scalar Fwtstar = `F_W_tstar'
		
end
